The relative benefits of environmental enrichment on learning and memory are greater when stressed: a meta-analysis of interactions in rodents
Supplementary Material S2
Setting-up
Loading packages
# devtools::install_github('Mikata-Project/ggthemr', force = TRUE)
pacman::p_load(tidyverse, here, metafor, clubSandwich, orchaRd, MuMIn, patchwork,
GoodmanKruskal, networkD3, ggplot2, ggsignif, visdat, ggalluvial, ggthemr, cowplot,
grDevices, png, grid, gridGraphics, pander, formatR)
# needed for model selection using MuMIn within metafor
eval(metafor:::.MuMIn)Loading data and functions
This loads the unprocessed data file and custom functions including
- calculating ‘focal’ effect sizes and variance for lnRR and SMD:
effect_set - calculating ‘pairwise’ comparisons (lnRR) and variance:
effect_set2 - calculating pair-wise contrasts between moderators:
contrast_fun
dat <- read_csv(here("Data", "Data_raw.csv"))
# Load custom function to extract data
source(here("R/Functions.R"))Data organisation
- removing study (Wang et al_2020) with negative values as lnRR cannot be calculated with negative values
- rounding down sample sizes that are reported as decimals due to averaging n across treatments
- getting effect sizes from function
- ‘flipping’ effect sizes so that all effect sizes are higher values = individuals do better in cognitive assays
- assigning human readable terms, and creating VCV of variance
# removing study with negative values as these are unable to be used for lnRR
dat <- droplevels(dat[!dat$First_author == 'Wang',])
#rounding down sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)
# 'Focal' effect_size
effect_size <- with(dat, mapply(effect_set,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size <- map_dfr(effect_size, I)
# 'Pairwise' effect size
effect_size2 <- with(dat, mapply(effect_set2,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size2 <- map_dfr(effect_size2, I)
full_info <- which(complete.cases(effect_size) == TRUE)
# adding effect sizes as column
dat <- bind_cols(dat, effect_size, effect_size2)
dat <- dat[full_info, ]
#Flipping 'lower is better' to 'higher is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E))
# currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))
# assigning human readable terms
dat <- dat %>% mutate(Type_assay = case_when(Type_assay == 1 ~ "Habituation",
Type_assay == 2 ~ "Conditioning",
Type_assay == 3 ~ "Recognition",
Type_assay == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
Learning_vs_memory == 2 ~ "Memory",
Learning_vs_memory == 3 ~ "Habituation"),
Type_reinforcement = case_when(Type_reinforcement== 1 ~"Appetitive",
Type_reinforcement== 2 ~ "Aversive",
Type_reinforcement== 3 ~ "Not applicable",
Type_reinforcement== 4 ~ "Unclear"),
Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
Type_stress_exposure == 2 ~ "Scent",
Type_stress_exposure == 3 ~ "Shock",
Type_stress_exposure == 4 ~ "Exertion",
Type_stress_exposure == 5 ~ "Restraint",
Type_stress_exposure == 6 ~ "MS",
Type_stress_exposure == 7 ~ "Circadian rhythm",
Type_stress_exposure == 8 ~ "Noise",
Type_stress_exposure == 9 ~ "Other",
Type_stress_exposure == 10 ~ "Combination",
Type_stress_exposure == 11 ~ "unclear"),
Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
Age_stress_exposure == 2 ~ "Early postnatal",
Age_stress_exposure == 3 ~ "Adolescent",
Age_stress_exposure == 4 ~ "Adult",
Age_stress_exposure == 5 ~ "Unclear"),
Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic",
Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"),
EE_social = case_when(EE_social == 1 ~ "Social",
EE_social== 2 ~ "Non-social",
EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise",
EE_exercise == 2 ~ "No exercise"),
Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Early postnatal",
Age_EE_exposure == 3 ~ "Adolescent",
Age_EE_exposure == 4 ~ "Adult",
Age_EE_exposure == 5 ~ "Unclear"),
Exposure_order = case_when(Exposure_order == 1 ~ "Stress first",
Exposure_order == 2 ~ "Enrichment first",
Exposure_order == 3 ~ "Concurrently",
Exposure_order == 4 ~ "Unclear"),
Age_assay = case_when(Age_assay == 1 ~ "Early postnatal",
Age_assay == 2 ~ "Adolescent",
Age_assay == 3 ~ "Adult",
Age_assay == 4 ~ "Unclear"),
Sex = case_when(Sex == 1 ~ "Female",
Sex == 2 ~ "Male",
Sex == 3 ~ "Mixed",
Sex == 4 ~ "Unclear"),
Type_EE_exposure = case_when(Type_EE_exposure == 1 ~ "Nesting material",
Type_EE_exposure == 2 ~ "Objects",
Type_EE_exposure == 3 ~ "Cage complexity",
Type_EE_exposure == 4 ~ "Wheel/trademill",
Type_EE_exposure == 5 ~ "Combination",
Type_EE_exposure == 6 ~ "Other",
Type_EE_exposure == 7 ~ "Unclear"),
ROB_blinding = case_when(ROB_blinding == 1 ~ "Yes",
ROB_blinding == 2 ~ "No",
ROB_blinding == 3 ~ "Unclear"),
ROB_randomisation = case_when(ROB_randomisation == 1 ~ "Yes",
ROB_randomisation == 2 ~ "No",
ROB_randomisation == 3 ~ "Unclear"))
#making variance VCVs
VCV_E <- impute_covariance_matrix(vi = dat$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
VCV_S <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ES <- impute_covariance_matrix(vi = dat$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
VCV_Ea <- impute_covariance_matrix(vi = dat$SMDV_E, cluster = dat$Study_ID, r = 0.5)
VCV_Sa <- impute_covariance_matrix(vi = dat$SMDV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ESa <- impute_covariance_matrix(vi = dat$SMDV_ES, cluster = dat$Study_ID, r = 0.5)
#write.csv(dat, file = here("Data", 'Data_processed.csv'), row.names = TRUE)Data exploration
General
Number of effect sizes:
# Number of effect sizes
length(unique(dat$ES_ID))## [1] 92
Number of studies:
length(unique(dat$Study_ID))## [1] 30
Publication year range:
min(dat$Year_published)## [1] 2006
max(dat$Year_published)## [1] 2021
Visual of missing data
plot_missing <- vis_miss(dat) + theme(plot.title = element_text(hjust = 0.5, vjust = 3),
plot.margin = margin(t = 0.5, r = 3, b = 1, l = 1, unit = "cm")) + ggtitle("Missing data in the selected predictors") #no missing values
plot_missing# useGoodman and Kruskal’s τ measure of association between categorical
# predictor variables (function from package GoodmanKruskal:
# https://cran.r-project.org/web/packages/GoodmanKruskal/vignettes/GoodmanKruskal.html)
# GKmatrix <- GKtauDataframe(subset(dat, select = c('Sex', 'Type_assay',
# 'Learning_vs_memory', #'Type_reinforcement', 'Type_stress_exposure',
# 'Age_stress_exposure', 'Stress_duration', #'EE_social_HR', 'EE_exercise',
# 'Age_EE_exposure', 'Exposure_order', 'Age_assay'))) plot(GKmatrix)
# simple pairwise contingency tables table(dat$Type_assay,
# dat$Type_reinforcement) table(dat$Age_stress_exposure, dat$Age_EE_exposure)
# table(dat$Type_stress_exposure, dat$Age_stress_exposure)
# table(dat$Type_stress_exposure, dat$Age_assay)
# table(dat$Type_stress_exposure, dat$Stress_duration)Visual of missing data for each column. Note that most missing data are for different descriptive statistics as papers are often report results using different types of descriptive statistics (i.e., SD, SE, mean, median etc). Cannot reduce text size in this figure
Alluvial diagrams
Shows the relationship/nestedness of different data elements (used to produce Fig. 2)
Subjects info: species-strain-sex
freq_A <- as.data.frame(table(dat$Sex, dat$Common_species, dat$Strain)) %>%
rename(Sex = Var1, Species = Var2, Strain = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_A), axes = 1:3, silent = TRUE)
p1 <- ggplot(data = freq_A, aes(axis1 = Sex, axis2 = Species, axis3 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Sex)) + geom_flow() + geom_stratum(aes(fill = Sex)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Sex",
"Species", "Strain"), expand = c(0.15, 0.05), position = "top") + ggtitle("A study subjects")
p1 Fig. 2a Relationship/nestedness between study species, strain of rodent, and sex
Stress info: age-duration-type stress
freq_C <- as.data.frame(table(dat$Age_stress_exposure, dat$Stress_duration, dat$Type_stress_exposure)) %>%
rename(Age_stress = Var1, Duration_stress = Var2, Type_stress = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_C), axes = 1:3, silent = TRUE)
p3 <- ggplot(data = freq_C, aes(axis1 = Age_stress, axis2 = Duration_stress, axis3 = Type_stress,
y = Freq)) + geom_alluvium(aes(fill = Age_stress)) + geom_flow() + geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Age",
"Duration", "Type"), expand = c(0.1, 0.1), position = "top") + ggtitle("C stress exposure")
p3 Fig. 2c Relationship/nestedess between the age of stress exposure, the duration of the stres, and the type of stress
Assay info: LvsM-type-reinforcement
freq_D <- as.data.frame(table(dat$Learning_vs_memory, dat$Type_assay, dat$Type_reinforcement)) %>%
rename(Learning_Memory = Var1, Type = Var2, Reinforcement = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_D), axes = 1:3, silent = TRUE)
p4 <- ggplot(data = freq_D, aes(axis1 = Learning_Memory, axis2 = Type, axis3 = Reinforcement,
y = Freq)) + geom_alluvium(aes(fill = Learning_Memory)) + geom_flow() + geom_stratum(aes(fill = Learning_Memory)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Learning_Memory",
"Type", "Reinforcement"), expand = c(0.1, 0.1), position = "top") + ggtitle("D cognitive assay")
p4 Fig. 2d Relationship/nestedness between if the assay measured a learning or memory trait, the type of assay used, and the type of reinforcement used
Combined plot
# p1 + scale_fill_brewer(palette = 'Set3') #Pastel1
(p1 + scale_fill_brewer(palette = "Set3"))/(p2 + scale_fill_brewer(palette = "Set3"))/(p3 +
scale_fill_brewer(palette = "Set3"))/(p4 + scale_fill_brewer(palette = "Set3")) +
plot_layout(ncol = 1, heights = c(1, 1, 1, 1, 1))# ggsave(file = './figs/Alluvial_diagrams_v2.pdf', width = 10, height = 12,
# units = 'cm', dpi = 300, scale = 2) #, device = cairo_pdf)Fig. 2 The combined plot of all alluvial diagrams used to make Fig. 2 in the main text
Risk of Bias
Percent of studies that used randomisation:
dat %>%
group_by(ROB_randomisation) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_randomisation n
## <chr> <int>
## 1 Unclear 15
## 2 Yes 15
15/30## [1] 0.5
Percent of studies that used blinding:
# blinding
dat %>%
group_by(ROB_blinding) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_blinding n
## <chr> <int>
## 1 Unclear 24
## 2 Yes 6
6/30## [1] 0.2
Modelling with lnRR
Main effect: environmental enrichment
Meta-analysis
To quantify differences in individual performance in cognitive assays with environmental enrichment, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{EE}} = \ln\left( \frac{ M_\text{ES} + M_\text{EC}} {2} \right) - \ln\left( \frac {M_\text{CS} + M_\text{CC}} {2} \right) = \ln \left( {\frac {M_\text{ES} + M_\text{EC}} {M_\text{CS} + M_\text{CC}}} \right), \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{EE}}) = \left( \frac{1} { M_\text{ES} + M_\text{EC} } \right)^2 \left( \frac{SD_\text{ES}^2}{N_\text{ES}} + \frac{SD_\text{EC}^2}{N_\text{EC}} \right) + \left( \frac {1} { M_\text{CS} + M_\text{CC} } \right)^2 \left(\frac{SD_\text{CS}^2}{N_\text{CS}} + \frac{SD_\text{CC}^2}{n_\text{CC}} \right) \] Variance was converted in to variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variance from the same studies.
#dat <- read_csv(here("Data","Data_processed.csv"))
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.7836 19.5672 27.5672 37.6106 28.0323
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0082 0.0905 30 no Study_ID
## sigma^2.2 0.0345 0.1858 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 809.2712, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1860 0.0320 5.8116 91 <.0001 0.1224 0.2496 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.310673e-01 1.786605e-01 7.524068e-01 6.097041e-10
orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the main effect of environmental enrichment model, we conducted a series of uni-moderator analyses. We calculated marginal R2 for each moderator as well as conducted a series of pair-wise contrasts between in moderator category.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the response learning or memory?
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.3793 18.7586 28.7586 41.2577 29.4729
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0076 0.0873 30 no Study_ID
## sigma^2.2 0.0340 0.1843 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2303 0.0450 5.1227 90 <.0001 0.1410
## Learning_vs_memoryMemory 0.1624 0.0355 4.5713 90 <.0001 0.0919
## ci.ub
## Learning_vs_memoryLearning 0.3197 ***
## Learning_vs_memoryMemory 0.2330 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E1) ## R2_marginal R2_coditional
## 0.02572583 1.00000000
LvsM_E<- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21) + # mean estimate
scale_size_continuous(range = c(1, 7)) + # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_E Fig. 4a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_E1 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = Learning_vs_memory,
VCV = VCV_E)
res_table_mod_E1 <- get_estimate(model = mod_E1, contra = contra_mod_E1, moderator = Learning_vs_memory)
res_table_mod_E1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning 0.230 0.141 0.320 0.00000170 -0.185 0.645
## 2 Memory 0.162 0.0919 0.233 0.0000154 -0.249 0.574
## 3 Learning-Memory -0.0679 -0.164 0.0284 0.165 -0.484 0.349
Type of assay
The broad category of the type of assay used to measure learning or memory
dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_E2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8496 15.6993 27.6993 42.6311 28.7237
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0995 30 no Study_ID
## sigma^2.2 0.0321 0.1792 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.2167 0.0351 6.1650 89 <.0001 0.1468 0.2865
## Type_assayHabituation 0.1821 0.1147 1.5886 89 0.1157 -0.0457 0.4100
## Type_assayRecognition 0.0554 0.0640 0.8659 89 0.3889 -0.0717 0.1826
##
## Type_assayConditioning ***
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E2) ## R2_marginal R2_coditional
## 0.07376134 1.00000000
Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Conditioning", "Recognition")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_E Fig. 4b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E2 <- contrast_fun(data = dat1, response = lnRR_Ea, moderator = Type_assay,
VCV = VCV_E1)
res_table_mod_E2 <- get_estimate(model = mod_E2, contra = contra_mod_E2, moderator = Type_assay)
res_table_mod_E2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning 0.217 0.147 0.287 2.02e-8 -0.197 0.630
## 2 Habituation 0.182 -0.0457 0.410 1.16e-1 -0.285 0.649
## 3 Recognition 0.0554 -0.0717 0.183 3.89e-1 -0.371 0.482
## 4 Conditioning-Habitua… -0.0345 -0.262 0.193 7.63e-1 -0.501 0.432
## 5 Conditioning-Recogni… -0.161 -0.296 -0.0265 1.96e-2 -0.590 0.268
## 6 Habituation-Recognit… -0.127 -0.381 0.127 3.24e-1 -0.607 0.353
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8702 15.7405 27.7405 42.6723 28.7649
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0109 0.1046 30 no Study_ID
## sigma^2.2 0.0319 0.1787 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.2175 0.0726 2.9942 89 0.0036 0.0732
## Type_reinforcementAversive 0.2191 0.0410 5.3456 89 <.0001 0.1377
## Type_reinforcementNot applicable 0.0812 0.0560 1.4504 89 0.1505 -0.0300
## ci.ub
## Type_reinforcementAppetitive 0.3618 **
## Type_reinforcementAversive 0.3006 ***
## Type_reinforcementNot applicable 0.1924
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E3) ## R2_marginal R2_coditional
## 0.07720103 1.00000000
Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Aversive", "Not applicable")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_E Fig. 4c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E3 <- contrast_fun(data = dat2, response = lnRR_Ea, moderator = Type_reinforcement,
VCV = VCV_E2)
res_table_mod_E3 <- get_estimate(model = mod_E3, contra = contra_mod_E3, moderator = Type_reinforcement)
res_table_mod_E3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive 0.217 0.0732 0.362 3.56e-3 -0.218 0.653
## 2 Aversive 0.219 0.138 0.301 6.89e-7 -0.200 0.638
## 3 Not applicable 0.0812 -0.0300 0.192 1.50e-1 -0.345 0.507
## 4 Appetitive-Aversive 0.00164 -0.163 0.167 9.84e-1 -0.442 0.445
## 5 Appetitive-Not applic… -0.136 -0.315 0.0425 1.33e-1 -0.585 0.312
## 6 Aversive-Not applicab… -0.138 -0.259 -0.0172 2.56e-2 -0.567 0.291
Age of enrichment
The age when individuals were exposed to environmental enrichment
dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)
mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_E4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.9490 13.8980 23.8980 36.1698 24.6480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0078 0.0885 29 no Study_ID
## sigma^2.2 0.0327 0.1808 88 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdolescent 0.1799 0.0370 4.8553 86 <.0001 0.1062 0.2535
## Age_EE_exposureAdult 0.2340 0.0620 3.7734 86 0.0003 0.1107 0.3573
##
## Age_EE_exposureAdolescent ***
## Age_EE_exposureAdult ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E4) ## R2_marginal R2_coditional
## 0.01127347 1.00000000
Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_E Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E4 <- contrast_fun(data = dat3, response = lnRR_Ea, moderator = Age_EE_exposure,
VCV = VCV_E3)
res_table_mod_E4 <- get_estimate(model = mod_E4, contra = contra_mod_E4, moderator = Age_EE_exposure)
res_table_mod_E4## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.180 0.106 0.254 0.00000533 -0.227 0.587
## 2 Adult 0.234 0.111 0.357 0.000295 -0.185 0.653
## 3 Adolescent-Adult 0.0541 -0.0895 0.198 0.456 -0.371 0.479
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.0303 20.0606 30.0606 42.5597 30.7749
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0096 0.0980 30 no Study_ID
## sigma^2.2 0.0345 0.1857 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1849 0.0407 4.5464 90 <.0001 0.1041 0.2657
## EE_exerciseNo exercise 0.1900 0.0556 3.4151 90 0.0010 0.0795 0.3005
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E5) ## R2_marginal R2_coditional
## 0.0001360993 0.9999999987
Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Exercise_E Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E5 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = EE_exercise,
VCV = VCV_E)
res_table_mod_E5 <- get_estimate(model = mod_E5, contra = contra_mod_E5, moderator = EE_exercise)
res_table_mod_E5## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Exercise 0.185 0.104 0.266 0.0000169 -0.240 0.610
## 2 No exercise 0.190 0.0795 0.301 0.000958 -0.241 0.621
## 3 Exercise-No exercise 0.00511 -0.132 0.142 0.941 -0.434 0.444
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"))
# VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster =
# dat_Efm$Study_ID, r = 0.5)
mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Efm)
# summary(mod_Efm) r2_ml(mod_Efm)
res_Efm <- dredge(mod_Efm, trace = 2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2 <- subset(res_Efm, delta <= 6, recalc.weights = FALSE)
importance(res_Efm2)## Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights: 0.63 0.36 0.16 0.14
## N containing models: 11 10 5 5
## Learning_vs_memory Type_reinforcement
## Sum of weights: 0.10 0.07
## N containing models: 4 4
Funnel plot
We produced funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry
Used to produce Fig. 7a
# funnel plot
Funnel_E <- funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")# Funnel_EFig.7a Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year
VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID,
r = 0.5)
dat_Efm$c_Year_published <- as.vector(scale(dat_Efm$Year_published, scale = F))
dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_E <- rma.mv(lnRR_Sa, VCV_Efm, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory +
Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise +
Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
test = "t", data = dat_Efm)
# PB_MR_E
estimates_PB_MR_E <- estimates.CI(PB_MR_E)
estimates_PB_MR_E## estimate mean lower upper
## 1 intrcpt -9.48980951 -50.89864510 31.91902609
## 2 sqrt_inv_e_n -0.06281057 -0.86857758 0.74295645
## 3 Learning_vs_memoryMemory -0.06858190 -0.17753633 0.04037252
## 4 Year_published 0.00467568 -0.01591232 0.02526368
## 5 Type_assayHabituation -0.55467648 -0.99644702 -0.11290594
## 6 Type_assayRecognition -0.11960487 -0.45315282 0.21394309
## 7 Type_reinforcementAversive 0.15536918 -0.15356143 0.46429979
## 8 Type_reinforcementNot applicable 0.32741017 -0.11084706 0.76566739
## 9 EE_socialSocial 0.04883161 -0.16642224 0.26408546
## 10 EE_exerciseNo exercise -0.03983196 -0.28097491 0.20131100
## 11 Age_stress_exposureAdult -0.16033122 -0.54043115 0.21976871
## 12 Age_stress_exposureEarly postnatal -0.05978986 -0.45195633 0.33237661
## 13 Age_stress_exposurePrenatal -0.13988902 -0.53943088 0.25965285
Leave-one-out analysis
We individually removed studies to determine if any singular studies hve a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Eb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = VCV_Eb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# TODO Shinichi to check
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0) {
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_E, file = here("Rdata", "MA_CVR_E.rds"))Used to produce Fig. S3
# telling ggplot to stop reordering factors
MA_CVR_E <- readRDS(file = here("Rdata", "MA_CVR_E.rds"))
MA_CVR_E$left_out <- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out <- factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)
# plotting
leaveoneout_E <- ggplot(MA_CVR_E) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_E0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Edat$Study_ID <- as.integer(dat$Study_ID)Fig. S3 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the main effect of environmental enrichment meta-analytic model using standardised mean difference (SMD) instead of lnRR
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -111.6284 223.2568 231.2568 241.3003 231.7220
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2512 0.5012 30 no Study_ID
## sigma^2.2 0.4247 0.6517 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7241 0.1295 5.5895 91 <.0001 0.4668 0.9814 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09
Risk of Bias
Including if the study was randomised as a moderator and conducting a pair-wise contrast
mod_randomised_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_randomisation - 1,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_E)
r2_ml(mod_randomised_E)contra_mod_randomised_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_randomisation,
VCV = VCV_E)
res_table_mod_randomised_E <- get_estimate(model = mod_randomised_E, contra = contra_mod_randomised_E,
moderator = ROB_randomisation)
res_table_mod_randomised_E## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.140 0.0488 0.232 0.00304 -0.276 0.556
## 2 Yes 0.227 0.139 0.315 0.00000158 -0.188 0.642
## 3 Unclear-Yes 0.0868 -0.0352 0.209 0.161 -0.337 0.511
Including if the study used blinding and conducting a pair-wise contrast
mod_blinding_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_E)
r2_ml(mod_blinding_E)contra_mod_blinding_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_blinding,
VCV = VCV_E)
res_table_mod_blinding_E <- get_estimate(model = mod_blinding_E, contra = contra_mod_blinding_E,
moderator = ROB_blinding)
res_table_mod_blinding_E## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.198 0.122 0.274 0.00000131 -0.226 0.622
## 2 Yes 0.148 0.00811 0.288 0.0384 -0.292 0.588
## 3 Unclear-Yes -0.0501 -0.208 0.108 0.530 -0.496 0.396
Main effect: stress
Meta-analysis
To quantify differences in individual performance in cognitive assays with stress, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{stress}} = \ln\left( \frac{ M_\text{ES} + M_\text{CS}} {2} \right) - \ln\left( \frac {M_\text{EC} + M_\text{CC}} {2} \right) = \ln \left( {\frac {M_\text{ES} + M_\text{CS}} {M_\text{EC} + M_\text{CC}}} \right), \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{stress}}) = \left( \frac{1} { M_\text{ES} + M_\text{CS} } \right)^2 \left( \frac{SD_\text{ES}^2}{N_\text{ES}} + \frac{SD_\text{CS}^2}{N_\text{CS}} \right) + \left( \frac {1} { M_\text{EC} + M_\text{CC} } \right)^2 \left(\frac{SD_\text{EC}^2}{N_\text{EC}} + \frac{SD_\text{CC}^2}{N_\text{CC}} \right) \] Variance was converted in to variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variance from the same studies.
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.1156 28.2312 36.2312 46.2747 36.6964
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0993 30 no Study_ID
## sigma^2.2 0.0377 0.1941 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 946.9234, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1052 0.0337 -3.1217 91 0.0024 -0.1721 -0.0383 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.376124e-01 1.946168e-01 7.429955e-01 1.975455e-10
orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the main effect of stress model, we conducted a series of uni-moderator analyses. We calculated marginal R2 for each moderator as well as conducted a series of pair-wise contrasts between in moderator category.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the response learning or memory?
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.5281 29.0562 39.0562 51.5553 39.7705
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0094 0.0970 30 no Study_ID
## sigma^2.2 0.0388 0.1969 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.1211 0.0476 -2.5423 90 0.0127 -0.2157
## Learning_vs_memoryMemory -0.0974 0.0380 -2.5648 90 0.0120 -0.1728
## ci.ub
## Learning_vs_memoryLearning -0.0265 *
## Learning_vs_memoryMemory -0.0219 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S1) ## R2_marginal R2_coditional
## 0.002776034 1.000000000
LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_S Fig. 5a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S1 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Learning_vs_memory,
VCV = VCV_S)
res_table_mod_S1 <- get_estimate(model = mod_S1, contra = contra_mod_S1, moderator = Learning_vs_memory)
res_table_mod_S1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning -0.121 -0.216 -0.0265 0.0127 -0.567 0.325
## 2 Memory -0.0974 -0.173 -0.0219 0.0120 -0.540 0.345
## 3 Learning-Memory 0.0237 -0.0780 0.125 0.644 -0.424 0.472
Type of assay
The broad category of the type of assay used to measure learning or memory
dat$Type_assay<-as.factor(dat$Type_assay)
VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_S2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.8028 19.6055 31.6055 46.5374 32.6299
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0161 0.1271 30 no Study_ID
## sigma^2.2 0.0279 0.1671 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning -0.0981 0.0375 -2.6192 89 0.0104 -0.1725 -0.0237
## Type_assayHabituation -0.4615 0.1126 -4.0969 89 <.0001 -0.6853 -0.2377
## Type_assayRecognition -0.0534 0.0645 -0.8287 89 0.4095 -0.1816 0.0747
##
## Type_assayConditioning *
## Type_assayHabituation ***
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S2) ## R2_marginal R2_coditional
## 0.1853359 1.0000000
Learning_S <-orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Habituation", "Recognition")),
annotations = "***") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_S Fig. 5b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S2 <- contrast_fun(data = dat1, response = lnRR_Sa, moderator = Type_assay,
VCV = VCV_S1)
res_table_mod_S2 <- get_estimate(model = mod_S2, contra = contra_mod_S2, moderator = Type_assay)
res_table_mod_S2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning -0.0981 -0.173 -0.0237 1.04e-2 -0.522 0.326
## 2 Habituation -0.461 -0.685 -0.238 9.21e-5 -0.935 0.0119
## 3 Recognition -0.0534 -0.182 0.0747 4.10e-1 -0.490 0.383
## 4 Conditioning-Habituation -0.363 -0.584 -0.142 1.55e-3 -0.835 0.109
## 5 Conditioning-Recognition 0.0447 -0.0882 0.177 5.06e-1 -0.393 0.482
## 6 Habituation-Recognition 0.408 0.161 0.655 1.48e-3 -0.0768 0.893
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.8810 27.7621 39.7621 54.6939 40.7864
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0103 0.1016 30 no Study_ID
## sigma^2.2 0.0387 0.1966 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
##
## Model Results:
##
## estimate se tval df pval
## Type_reinforcementAppetitive -0.1846 0.0749 -2.4649 89 0.0156
## Type_reinforcementAversive -0.0730 0.0427 -1.7081 89 0.0911
## Type_reinforcementNot applicable -0.1172 0.0590 -1.9851 89 0.0502
## ci.lb ci.ub
## Type_reinforcementAppetitive -0.3334 -0.0358 *
## Type_reinforcementAversive -0.1579 0.0119 .
## Type_reinforcementNot applicable -0.2345 0.0001 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S3) ## R2_marginal R2_coditional
## 0.0366428 1.0000000
Reinforcement_S <-orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_S Fig. 5c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S3 <- contrast_fun(data = dat2, response = lnRR_Sa, moderator = Type_reinforcement,
VCV = VCV_S2)
res_table_mod_S3 <- get_estimate(model = mod_S3, contra = contra_mod_S3, moderator = Type_reinforcement)
res_table_mod_S3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive -0.185 -0.333 -3.58e-2 0.0156 -0.649 0.280
## 2 Aversive -0.0730 -0.158 1.19e-2 0.0911 -0.521 0.375
## 3 Not applicable -0.117 -0.235 1.12e-4 0.0502 -0.572 0.338
## 4 Appetitive-Aversive 0.112 -0.0592 2.82e-1 0.198 -0.360 0.583
## 5 Appetitive-Not applicable 0.0674 -0.119 2.54e-1 0.474 -0.410 0.545
## 6 Aversive-Not applicable -0.0442 -0.173 8.46e-2 0.497 -0.502 0.414
Age of stress
The age when individuals were exposed to stress
mod_S4 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.4083 24.8166 38.8166 56.1579 40.2166
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0048 0.0694 30 no Study_ID
## sigma^2.2 0.0392 0.1979 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0074 0.1159 0.0641 88 0.9490
## Age_stress_exposureAdult -0.2279 0.0622 -3.6664 88 0.0004
## Age_stress_exposureEarly postnatal -0.0561 0.0435 -1.2891 88 0.2007
## Age_stress_exposurePrenatal -0.1145 0.0743 -1.5404 88 0.1271
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.2228 0.2377
## Age_stress_exposureAdult -0.3514 -0.1044 ***
## Age_stress_exposureEarly postnatal -0.1425 0.0304
## Age_stress_exposurePrenatal -0.2621 0.0332
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4) ## R2_marginal R2_coditional
## 0.09987307 1.00000000
Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_signif(comparisons = list(c("Adult", "Early postnatal")),
annotations = "*") +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_S Fig. 5d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S4 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Age_stress_exposure,
VCV = VCV_S)
res_table_mod_S4 <- get_estimate(model = mod_S4, contra = contra_mod_S4, moderator = Age_stress_exposure)
res_table_mod_S4## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.00743 -0.223 0.238 9.49e-1 -0.469 0.484
## 2 Adult -0.228 -0.351 -0.104 4.20e-4 -0.663 0.207
## 3 Early postnatal -0.0561 -0.142 0.0304 2.01e-1 -0.482 0.370
## 4 Prenatal -0.114 -0.262 0.0332 1.27e-1 -0.557 0.328
## 5 Adolescent-Adult -0.235 -0.497 0.0260 7.70e-2 -0.727 0.257
## 6 Adolescent-Early postna… -0.0635 -0.309 0.182 6.09e-1 -0.547 0.421
## 7 Adolescent-Prenatal -0.122 -0.395 0.152 3.78e-1 -0.620 0.377
## 8 Adult-Early postnatal 0.172 0.0211 0.323 2.60e-2 -0.271 0.615
## 9 Adult-Prenatal 0.113 -0.0791 0.306 2.45e-1 -0.346 0.573
## 10 Early postnatal-Prenatal -0.0584 -0.229 0.113 4.99e-1 -0.509 0.392
Type of stress
The type of stressor used
dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)
mod_S5 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_S5) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.4138 22.8276 36.8276 53.5887 38.3618
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0115 0.1071 25 no Study_ID
## sigma^2.2 0.0401 0.2002 85 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination -0.0500 0.0892 -0.5605 81 0.5767 -0.2274
## Type_stress_exposureMS -0.0539 0.0560 -0.9630 81 0.3384 -0.1653
## Type_stress_exposureNoise -0.1203 0.1036 -1.1608 81 0.2491 -0.3265
## Type_stress_exposureRestraint -0.2101 0.0704 -2.9863 81 0.0037 -0.3501
## ci.ub
## Type_stress_exposureCombination 0.1274
## Type_stress_exposureMS 0.0575
## Type_stress_exposureNoise 0.0859
## Type_stress_exposureRestraint -0.0701 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S5)## R2_marginal R2_coditional
## 0.07029554 1.00000000
Stressor<- orchard_plot(mod_S5, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Stressor Fig. 5e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S5 <- contrast_fun(data = dat5, response = lnRR_Sa, moderator = Type_stress_exposure,
VCV = VCV_S3)
res_table_mod_S5 <- get_estimate(model = mod_S5, contra = contra_mod_S5, moderator = Type_stress_exposure)
res_table_mod_S5## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Combination -0.0500 -0.227 0.127 0.577 -0.535 0.435
## 2 MS -0.0539 -0.165 0.0575 0.338 -0.519 0.411
## 3 Noise -0.120 -0.326 0.0859 0.249 -0.617 0.376
## 4 Restraint -0.210 -0.350 -0.0701 0.00373 -0.683 0.263
## 5 Combination-MS -0.00393 -0.213 0.206 0.970 -0.502 0.494
## 6 Combination-Noise -0.0703 -0.342 0.202 0.608 -0.598 0.457
## 7 Combination-Restraint -0.160 -0.386 0.0658 0.162 -0.665 0.345
## 8 MS-Noise -0.0664 -0.301 0.168 0.575 -0.575 0.443
## 9 MS-Restraint -0.156 -0.335 0.0227 0.0861 -0.642 0.330
## 10 Noise-Restraint -0.0898 -0.339 0.159 0.475 -0.606 0.426
Stess duration
Was the stress acute or chronic?
dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)
mod_S6 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_S6) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.7898 27.5796 37.5796 49.9092 38.3204
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0062 0.0786 29 no Study_ID
## sigma^2.2 0.0391 0.1978 89 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.0135 0.0659 0.2045 87 0.8384 -0.1176 0.1445
## Stress_durationChronic -0.1360 0.0373 -3.6456 87 0.0005 -0.2101 -0.0618
##
## Stress_durationAcute
## Stress_durationChronic ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S6) ## R2_marginal R2_coditional
## 0.08245896 1.00000000
Duration_S <- orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
geom_signif(comparisons = list(c("Chronic", "Acute")),
annotations = "*") +
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Duration_S Fig. 5f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_S6 <- contrast_fun(data = dat6, response = lnRR_Sa, moderator = Stress_duration,
VCV = VCV_S4)
res_table_mod_S6 <- get_estimate(model = mod_S6, contra = contra_mod_S6, moderator = Stress_duration)
res_table_mod_S6## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Acute 0.0135 -0.118 0.145 0.838 -0.429 0.456
## 2 Chronic -0.136 -0.210 -0.0618 0.000453 -0.565 0.294
## 3 Acute-Chronic -0.149 -0.300 0.00111 0.0517 -0.599 0.300
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
# selecting moderator levels with k >=5
dat_Sfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), Type_stress_exposure %in%
c("Restraint", "Noise", "MS", "Combination"), Stress_duration %in% c("Chronic",
"Acute"))
# VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster =
# dat_Sfm$Study_ID, r = 0.5)
mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Sfm)
# summary(mod_Sfm) r2_ml(mod_Sfm)
res_Sfm <- dredge(mod_Sfm, trace = 2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2 <- subset(res_Sfm, delta <= 6, recalc.weights = FALSE)
importance(res_Sfm2)## Type_assay Stress_duration Type_stress_exposure
## Sum of weights: 0.91 0.88 0.22
## N containing models: 7 6 2
## Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights: 0.16 0.04 0.04
## N containing models: 2 1 1
Funnel plot
We produced funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.
Used to produce Fig. 7b
# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")# Funnel_SFig.7b Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.
# calculating inv effective sample size for use in full meta-regression
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID,
r = 0.5)
dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
# time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))
PB_MR_S <- rma.mv(lnRR_Sa, VCV_Sfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
Age_stress_exposure + Stress_duration, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), method = "REML", test = "t", data = dat_Sfm, control = list(optimizer = "optim",
optmethod = "Nelder-Mead"))
# PB_MR_S
estimates_PB_MR_S <- estimates.CI(PB_MR_S)
estimates_PB_MR_S## estimate mean lower upper
## 1 intrcpt 0.208757092 -0.47118467 0.88869885
## 2 sqrt_inv_e_n 0.109877719 -0.90636737 1.12612281
## 3 c_Year_published -0.006425821 -0.02804904 0.01519740
## 4 Type_assayHabituation -0.663428076 -1.09349872 -0.23335743
## 5 Type_assayRecognition -0.069789046 -0.41156145 0.27198335
## 6 Learning_vs_memoryMemory -0.076712961 -0.19562432 0.04219839
## 7 Type_reinforcementAversive 0.039563987 -0.13761334 0.21674132
## 8 Type_reinforcementNot applicable 0.189481532 -0.16955283 0.54851589
## 9 Type_stress_exposureMS 0.233030567 -0.40478886 0.87084999
## 10 Type_stress_exposureNoise -0.020700600 -0.76667312 0.72527192
## 11 Type_stress_exposureRestraint -0.067490308 -0.43366600 0.29868539
## 12 Age_stress_exposureAdult -0.129875123 -0.51040393 0.25065368
## 13 Age_stress_exposureEarly postnatal -0.193461086 -0.89857848 0.51165631
## 14 Age_stress_exposurePrenatal -0.146406300 -0.57377820 0.28096560
## 15 Stress_durationChronic -0.395711964 -0.65566557 -0.13575836
Leave-one-out sensitivity analysis
We individually removed studies to determine if any singular studies hve a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Sb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = VCV_Sb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_S0) {
df <- data.frame(est = mod_S0$b, lower = mod_S0$ci.lb, upper = mod_S0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_S, file = here("Rdata", "MA_CVR_S.rds"))Used to produce Fig. S4
MA_CVR_S <- readRDS(file = here("Rdata", "MA_CVR_S.rds"))
# telling ggplot to stop reordering factors
MA_CVR_S$left_out <- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out <- factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)
# plotting
leaveoneout_S <- ggplot(MA_CVR_S) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_S0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Sdat$Study_ID <- as.integer(dat$Study_ID)Fig. S4 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the main effect of stress meta-analytic model using standardised mean difference (SMD) instead of lnRR
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_S0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -120.2817 240.5634 248.5634 258.6068 249.0285
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3874 0.6224 30 no Study_ID
## sigma^2.2 0.4969 0.7049 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.4256 0.1498 -2.8403 91 0.0056 -0.7232 -0.1279 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09
Risk of Bias
Including if the study was randomised as a moderator and conducting a pair-wise contrast
mod_randomised_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_randomisation - 1,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_S)
r2_ml(mod_randomised_S)contra_mod_randomised_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_randomisation,
VCV = VCV_S)
res_table_mod_randomised_S <- get_estimate(model = mod_randomised_S, contra = contra_mod_randomised_S,
moderator = ROB_randomisation)
res_table_mod_randomised_S## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear -0.0672 -0.164 0.0291 0.169 -0.509 0.374
## 2 Yes -0.139 -0.230 -0.0480 0.00317 -0.580 0.301
## 3 Unclear-Yes -0.0719 -0.204 0.0607 0.284 -0.523 0.379
Including if the study used blinding and conducting a pair-wise contrast
mod_blinding_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_S)
r2_ml(mod_blinding_S)contra_mod_blinding_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_blinding,
VCV = VCV_S)
res_table_mod_blinding_S <- get_estimate(model = mod_blinding_S, contra = contra_mod_blinding_S,
moderator = ROB_blinding)
res_table_mod_blinding_S## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear -0.119 -0.196 -0.0419 0.00288 -0.564 0.326
## 2 Yes -0.0567 -0.202 0.0890 0.442 -0.518 0.405
## 3 Unclear-Yes 0.0625 -0.102 0.227 0.454 -0.406 0.531
Interaction: enrichment x stress
Meta-analysis
To quantify differences in individual performance in cognitive assay with the interaction of environmental enrichment and stress, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{interaction}} = \left(\ln M_\text{ES} - \ln M_\text{CS} \right) - \left( \ln M_\text{EC} - \ln M_\text{CC} \right) = \ln \left( \frac {M_\text{ES} / M_\text{CS} } {M_\text{EC} / \ M_\text{CC}} \right) \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{interaction}}) = \frac {SD_\text{ES}^2} { N_\text{ES} M_\text{ES}^2 } + \frac{SD_\text{EC}^2}{N_\text{EC}M_\text{EC}^2} + \frac {SD_\text{CS}^2}{N_\text{CS} M_\text{CS}^2} + \frac{SD_\text{CC}^2}{N_\text{CC} M_\text{CC}^2}. \] Variance was converted in to variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variance from the same studies.
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.8178 81.6355 89.6355 99.6790 90.1006
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0229 0.1513 92 no ES_ID
## sigma^2.3 0.0030 0.0544 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 303.2179, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1229 0.0596 2.0605 91 0.0422 0.0044 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.81703063 0.44913873 0.32576306 0.04212884
orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the interaction of environmental enrichment and stress, we conducted a series of uni-moderator analyses. We calculated marginal R2 for each moderator as well as conducted a series of pair-wise contrasts between in moderator category.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the response learning or memory?
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4769 80.9539 90.9539 103.4529 91.6682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0292 0.1708 30 no Study_ID
## sigma^2.2 0.0232 0.1524 92 no ES_ID
## sigma^2.3 0.0034 0.0582 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1744 0.0722 2.4166 90 0.0177 0.0310
## Learning_vs_memoryMemory 0.1057 0.0619 1.7065 90 0.0914 -0.0174
## ci.ub
## Learning_vs_memoryLearning 0.3178 *
## Learning_vs_memoryMemory 0.2288 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES1) ## R2_marginal R2_coditional
## 0.0197648 0.9405080
LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
LvsM_ES Fig. 6a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES1 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Learning_vs_memory,
VCV = VCV_ES)
res_table_mod_ES1 <- get_estimate(model = mod_ES1, contra = contra_mod_ES1, moderator = Learning_vs_memory)
res_table_mod_ES1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning 0.174 0.0310 0.318 0.0177 -0.316 0.665
## 2 Memory 0.106 -0.0174 0.229 0.0914 -0.379 0.591
## 3 Learning-Memory -0.0687 -0.176 0.0382 0.205 -0.550 0.413
Type of assay
The broad category of the type of assay used to measure learning or memory
VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_ES2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0874 78.1748 90.1748 105.1066 91.1992
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0370 0.1924 30 no Study_ID
## sigma^2.2 0.0192 0.1386 92 no ES_ID
## sigma^2.3 0.0018 0.0422 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.1525 0.0589 2.5877 89 0.0113 0.0354 0.2696
## Type_assayHabituation 0.1990 0.1415 1.4070 89 0.1629 -0.0820 0.4801
## Type_assayRecognition -0.0048 0.0800 -0.0606 89 0.9518 -0.1637 0.1541
##
## Type_assayConditioning *
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES2) ## R2_marginal R2_coditional
## 0.05775809 0.97105550
Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Conditioning", "Recognition")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Learning_ES Fig. 6b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES2 <- contrast_fun(data = dat1, response = lnRR_ESa, moderator = Type_assay,
VCV = VCV_ES1)
res_table_mod_ES2 <- get_estimate(model = mod_ES2, contra = contra_mod_ES2, moderator = Type_assay)
res_table_mod_ES2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning 0.153 0.0354 0.270 0.0113 -0.340 0.645
## 2 Habituation 0.199 -0.0820 0.480 0.163 -0.356 0.754
## 3 Recognition -0.00485 -0.164 0.154 0.952 -0.509 0.499
## 4 Conditioning-Habituation 0.0465 -0.219 0.312 0.728 -0.501 0.594
## 5 Conditioning-Recognition -0.157 -0.304 -0.0105 0.0360 -0.658 0.343
## 6 Habituation-Recognition -0.204 -0.488 0.0801 0.157 -0.760 0.353
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID, r = 0.5)
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0604 78.1208 90.1208 105.0526 91.1452
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0382 0.1954 30 no Study_ID
## sigma^2.2 0.0189 0.1377 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1007 0.1075 0.9366 89 0.3515 -0.1129
## Type_reinforcementAversive 0.1573 0.0569 2.7618 89 0.0070 0.0441
## Type_reinforcementNot applicable 0.0147 0.0702 0.2101 89 0.8341 -0.1247
## ci.ub
## Type_reinforcementAppetitive 0.3143
## Type_reinforcementAversive 0.2704 **
## Type_reinforcementNot applicable 0.1541
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES3) ## R2_marginal R2_coditional
## 0.0586952 1.0000000
Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Aversive", "Not applicable")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Reinforcement_ES Fig. 6c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES3 <- contrast_fun(data = dat2, response = lnRR_ESa, moderator = Type_reinforcement,
VCV = VCV_ES2)
res_table_mod_ES3 <- get_estimate(model = mod_ES3, contra = contra_mod_ES3, moderator = Type_reinforcement)
res_table_mod_ES3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive 0.101 -0.113 0.314 0.351 -0.420 0.621
## 2 Aversive 0.157 0.0441 0.270 0.00698 -0.331 0.645
## 3 Not applicable 0.0147 -0.125 0.154 0.834 -0.480 0.510
## 4 Appetitive-Aversive 0.0566 -0.182 0.295 0.639 -0.475 0.588
## 5 Appetitive-Not applicable -0.0859 -0.333 0.161 0.492 -0.621 0.450
## 6 Aversive-Not applicable -0.143 -0.278 -0.00744 0.0389 -0.636 0.351
Age of enrichment
The age when individuals were exposed to environmental enrichment
VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID, r = 0.5)
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_ES4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -36.6407 73.2813 83.2813 95.5531 84.0313
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0332 0.1822 29 no Study_ID
## sigma^2.2 0.0207 0.1439 88 no ES_ID
## sigma^2.3 0.0007 0.0265 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.1284 0.0583 2.2006 86 0.0304 0.0124
## Age_EE_exposureAdult 0.1247 0.0921 1.3538 86 0.1793 -0.0584
## ci.ub
## Age_EE_exposureAdolescent 0.2444 *
## Age_EE_exposureAdult 0.3077
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES4) ## R2_marginal R2_coditional
## 0.0000400928 0.9871095822
Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_enrichment_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES4 <- contrast_fun(data = dat3, response = lnRR_ESa, moderator = Age_EE_exposure,
VCV = VCV_ES3)
res_table_mod_ES4 <- get_estimate(model = mod_ES4, contra = contra_mod_ES4, moderator = Age_EE_exposure)
res_table_mod_ES4## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.128 0.0124 0.244 0.0304 -0.350 0.607
## 2 Adult 0.125 -0.0584 0.308 0.179 -0.375 0.624
## 3 Adolescent-Adult -0.00373 -0.213 0.205 0.972 -0.513 0.506
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_ES5<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4968 80.9935 90.9935 103.4926 91.7078
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0325 0.1803 30 no Study_ID
## sigma^2.2 0.0230 0.1517 92 no ES_ID
## sigma^2.3 0.0065 0.0805 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 297.3270, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.8401, p-val = 0.1647
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1051 0.0780 1.3474 90 0.1812 -0.0499 0.2602
## EE_exerciseNo exercise 0.1687 0.0967 1.7448 90 0.0844 -0.0234 0.3609
##
## EE_exerciseExercise
## EE_exerciseNo exercise .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES5) ## R2_marginal R2_coditional
## 0.01474459 0.89712469
Exercise_ES <- orchard_plot(mod_ES5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Exercise_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES5 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = EE_exercise,
VCV = VCV_ES)
res_table_mod_ES5 <- get_estimate(model = mod_ES5, contra = contra_mod_ES5, moderator = EE_exercise)
res_table_mod_ES5## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Exercise 0.105 -0.0499 0.260 0.181 -0.413 0.624
## 2 No exercise 0.169 -0.0234 0.361 0.0844 -0.362 0.699
## 3 Exercise-No exercise 0.0636 -0.138 0.265 0.532 -0.470 0.598
Age of stress
The age when individuals were exposed to stress
mod_ES7 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES7) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0917 78.1834 92.1834 109.5247 93.5834
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0295 0.1718 30 no Study_ID
## sigma^2.2 0.0232 0.1522 92 no ES_ID
## sigma^2.3 0.0091 0.0954 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 286.4225, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.5932, p-val = 0.1832
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent -0.0137 0.1762 -0.0775 88 0.9384
## Age_stress_exposureAdult 0.1677 0.1140 1.4708 88 0.1449
## Age_stress_exposureEarly postnatal 0.1067 0.0920 1.1602 88 0.2491
## Age_stress_exposurePrenatal 0.3179 0.1311 2.4254 88 0.0173
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.3639 0.3366
## Age_stress_exposureAdult -0.0589 0.3942
## Age_stress_exposureEarly postnatal -0.0761 0.2896
## Age_stress_exposurePrenatal 0.0574 0.5784 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES7) ## R2_marginal R2_coditional
## 0.09276574 0.86630830
Age_stress_ES<-orchard_plot(mod_ES7, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_stress_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES7 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Age_stress_exposure,
VCV = VCV_ES)
res_table_mod_ES7 <- get_estimate(model = mod_ES7, contra = contra_mod_ES7, moderator = Age_stress_exposure)
res_table_mod_ES7## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent -0.0137 -0.364 0.337 0.938 -0.619 0.592
## 2 Adult 0.168 -0.0589 0.394 0.145 -0.376 0.711
## 3 Early postnatal 0.107 -0.0761 0.290 0.249 -0.420 0.634
## 4 Prenatal 0.318 0.0574 0.578 0.0173 -0.241 0.876
## 5 Adolescent-Adult 0.181 -0.236 0.598 0.390 -0.465 0.828
## 6 Adolescent-Early postna… 0.120 -0.275 0.515 0.546 -0.512 0.753
## 7 Adolescent-Prenatal 0.332 -0.105 0.768 0.135 -0.328 0.991
## 8 Adult-Early postnatal -0.0609 -0.280 0.158 0.581 -0.601 0.479
## 9 Adult-Prenatal 0.150 -0.132 0.432 0.292 -0.419 0.719
## 10 Early postnatal-Prenatal 0.211 -0.0462 0.469 0.107 -0.346 0.768
Type of stress
The type of stressor used
VCV_ES5 <- impute_covariance_matrix(vi = dat5$lnRRV_ES, cluster = dat5$Study_ID, r = 0.5)
mod_ES8 <- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -34.4046 68.8091 82.8091 99.5703 84.3434
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0426 0.2064 25 no Study_ID
## sigma^2.2 0.0232 0.1524 85 no ES_ID
## sigma^2.3 0.0104 0.1021 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 281.9708, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 0.5137, p-val = 0.7258
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination 0.1111 0.1458 0.7618 81 0.4484 -0.1790
## Type_stress_exposureMS 0.1185 0.1099 1.0784 81 0.2840 -0.1001
## Type_stress_exposureNoise 0.1651 0.1795 0.9198 81 0.3604 -0.1920
## Type_stress_exposureRestraint 0.1374 0.1252 1.0978 81 0.2755 -0.1116
## ci.ub
## Type_stress_exposureCombination 0.4011
## Type_stress_exposureMS 0.3370
## Type_stress_exposureNoise 0.5221
## Type_stress_exposureRestraint 0.3865
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES8)## R2_marginal R2_coditional
## 0.004455703 0.863910284
Stressor_ES <- orchard_plot(mod_ES8, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Stressor_ES Fig. 6e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES8 <- contrast_fun(data = dat5, response = lnRR_ESa, moderator = Type_stress_exposure,
VCV = VCV_ES5)
res_table_mod_ES8 <- get_estimate(model = mod_ES8, contra = contra_mod_ES8, moderator = Type_stress_exposure)
res_table_mod_ES8## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Combination 0.111 -0.179 0.401 0.448 -0.510 0.732
## 2 MS 0.118 -0.100 0.337 0.284 -0.473 0.710
## 3 Noise 0.165 -0.192 0.522 0.360 -0.490 0.820
## 4 Restraint 0.137 -0.112 0.386 0.276 -0.466 0.741
## 5 Combination-MS 0.00741 -0.299 0.313 0.962 -0.621 0.636
## 6 Combination-Noise 0.0540 -0.352 0.460 0.792 -0.629 0.737
## 7 Combination-Restraint 0.0263 -0.300 0.353 0.873 -0.613 0.665
## 8 MS-Noise 0.0466 -0.312 0.405 0.797 -0.609 0.703
## 9 MS-Restraint 0.0189 -0.246 0.284 0.887 -0.591 0.629
## 10 Noise-Restraint -0.0276 -0.403 0.348 0.884 -0.693 0.638
Stress duration
Was the stress acute or chronic?
VCV_ES6 <- impute_covariance_matrix(vi = dat6$lnRRV_ES, cluster = dat6$Study_ID, r = 0.5)
mod_ES9 <-rma.mv(yi = lnRR_ESa, V = VCV_ES6, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_ES9) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -35.9010 71.8020 81.8020 94.1315 82.5427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0133 0.1155 29 no Study_ID
## sigma^2.2 0.0260 0.1611 89 no ES_ID
## sigma^2.3 0.0080 0.0895 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 278.4877, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.3877, p-val = 0.0153
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.0367 0.0930 -0.3946 87 0.6941 -0.2216 0.1482
## Stress_durationChronic 0.1854 0.0732 2.5347 87 0.0130 0.0400 0.3308
##
## Stress_durationAcute
## Stress_durationChronic *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES9) ## R2_marginal R2_coditional
## 0.1598311 0.8575918
Duration_ES<- orchard_plot(mod_ES9, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_signif(comparisons = list(c("Acute", "Chronic")),
annotations = "*") +
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Duration_ES Fig. 6f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES9 <- contrast_fun(data = dat6, response = lnRR_ESa, moderator = Stress_duration,
VCV = VCV_ES6)
res_table_mod_ES9 <- get_estimate(model = mod_ES9, contra = contra_mod_ES9, moderator = Stress_duration)
res_table_mod_ES9## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Acute -0.0367 -0.222 0.148 0.694 -0.507 0.433
## 2 Chronic 0.185 0.0400 0.331 0.0130 -0.271 0.642
## 3 Acute-Chronic 0.222 0.0381 0.406 0.0186 -0.248 0.692
Order to treatment exposure
The order in which individuals were exposed to enrichment and stress
mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order -1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES10)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0408 78.0817 90.0817 105.0135 91.1061
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0227 0.1507 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 292.2561, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1546, p-val = 0.0287
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Exposure_orderConcurrently 0.1492 0.1208 1.2351 89 0.2201 -0.0909
## Exposure_orderEnrichment first -0.1782 0.1659 -1.0744 89 0.2856 -0.5079
## Exposure_orderStress first 0.1370 0.0526 2.6046 89 0.0108 0.0325
## ci.ub
## Exposure_orderConcurrently 0.3893
## Exposure_orderEnrichment first 0.1514
## Exposure_orderStress first 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES10)## R2_marginal R2_coditional
## 0.1027207 1.0000000
Order_ES <- orchard_plot(mod_ES10, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Order_ES Fig. 6g Orchard plot showing the group-wise means of the categorical variable ‘Exposure_order’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_ES10 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Exposure_order,
VCV = VCV_ES)
res_table_mod_ES10 <- get_estimate(model = mod_ES10, contra = contra_mod_ES10, moderator = Exposure_order)
res_table_mod_ES10## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Concurrently 0.149 -0.0909 0.389 0.220 -0.372 0.671
## 2 Enrichment first -0.178 -0.508 0.151 0.286 -0.747 0.390
## 3 Stress first 0.137 0.0325 0.241 0.0108 -0.338 0.612
## 4 Concurrently-Enrichment … -0.327 -0.735 0.0803 0.114 -0.944 0.289
## 5 Concurrently-Stress first -0.0123 -0.274 0.250 0.926 -0.544 0.520
## 6 Enrichment first-Stress … 0.315 -0.0306 0.661 0.0735 -0.263 0.893
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
dat_ESfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"), Type_stress_exposure %in%
c("Restraint", "Noise", "MS", "Combination"), Stress_duration %in% c("Chronic",
"Acute"), Age_stress_exposure %in% c("Prenatal", "Early postnatal", "Adult"))
# VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster =
# dat_ESfm$Study_ID, r = 0.5)
mod_ESfm <- rma.mv(yi = lnRR_Sa, V = VCV_ESfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure + Type_stress_exposure +
Age_stress_exposure + Stress_duration + Exposure_order, random = list(~1 | Study_ID,
~1 | ES_ID, ~1 | Strain), test = "t", data = dat_ESfm)
# summary(mod_ESfm) r2_ml(mod_ESfm)
res_ESfm <- dredge(mod_ESfm, trace = 2)
saveRDS(res_ESfm, file = here("Rdata", "res_ESfm.rds"))
# also saving the full model and data
saveRDS(mod_ESfm, file = here("Rdata", "mod_ESfm.rds"))
saveRDS(dat_ESfm, file = here("Rdata", "dat_ESfm.rds"))The akaike weights for the top set of models with AIC < 6
dat_ESfm <- readRDS(file = here("Rdata", "dat_ESfm.rds"))
mod_ESfm <- readRDS(file = here("Rdata", "mod_ESfm.rds"))
res_ESfm <- readRDS(file = here("Rdata", "res_ESfm.rds"))
res_ESfm2 <- subset(res_ESfm, delta <= 6, recalc.weights = FALSE)
importance(res_ESfm2)## Type_assay Stress_duration Age_EE_exposure EE_exercise
## Sum of weights: 0.65 0.63 0.26 0.11
## N containing models: 19 18 7 6
## Learning_vs_memory EE_social Age_stress_exposure
## Sum of weights: 0.08 0.08 0.05
## N containing models: 4 4 2
## Type_stress_exposure Exposure_order
## Sum of weights: 0.03 0.02
## N containing models: 2 1
Funnel plot
We produced funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.
Used to produce Fig. 7c
Funnel_ES <- funnel(mod_ESfm, xlab = "lnRR", ylab = "Standard Error")# Funnel_ESFig.7c Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.
# calculating inv effective sample size for use in full meta-regression
dat_ESfm$sqrt_inv_e_n <- with(dat_ESfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
dat_ESfm$c_Year_published <- as.vector(scale(dat_ESfm$Year_published, scale = F))
VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster = dat_ESfm$Study_ID,
r = 0.5)
# time lag bias and eggers regression
PB_MR_ES <- rma.mv(lnRR_ESa, VCV_ESfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
Age_stress_exposure + EE_social + EE_exercise + Stress_duration + Exposure_order,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", test = "t",
data = dat_ESfm, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
estimates_PB_MR_ES <- estimates.CI(PB_MR_ES)
estimates_PB_MR_ES## estimate mean lower upper
## 1 intrcpt 0.104890890 -1.32422406 1.53400584
## 2 sqrt_inv_e_n -0.380617946 -2.34738544 1.58614955
## 3 c_Year_published -0.003614793 -0.03254341 0.02531382
## 4 Type_assayHabituation 0.225637592 -0.28343491 0.73471009
## 5 Type_assayRecognition 0.014042786 -0.37390986 0.40199543
## 6 Learning_vs_memoryMemory 0.027907733 -0.10104526 0.15686073
## 7 Type_reinforcementAversive 0.074150779 -0.20556264 0.35386420
## 8 Type_reinforcementNot applicable -0.180164458 -0.64972276 0.28939384
## 9 Type_stress_exposureMS -0.670185502 -1.72077062 0.38039961
## 10 Type_stress_exposureNoise -0.307675466 -1.41676992 0.80141899
## 11 Type_stress_exposureRestraint -0.159808179 -0.87246493 0.55284857
## 12 Age_stress_exposureEarly postnatal 0.472151143 -0.39118506 1.33548735
## 13 Age_stress_exposurePrenatal 0.199706727 -0.39069154 0.79010499
## 14 EE_socialSocial 0.134840904 -0.25767695 0.52735875
## 15 EE_exerciseNo exercise 0.127498329 -0.18382738 0.43882404
## 16 Stress_durationChronic 0.449057668 0.05381018 0.84430516
## 17 Exposure_orderEnrichment first -0.132255288 -0.74038190 0.47587133
## 18 Exposure_orderStress first -0.057436396 -0.54737789 0.43250509
Leave-one-out analysis
We individually removed studies to determine if any singular studies hve a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_ESb <- impute_covariance_matrix(vi = d$lnRRV_ES, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_ESa, V = VCV_ESb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_ES0) {
df <- data.frame(est = mod_ES0$b, lower = mod_ES0$ci.lb, upper = mod_ES0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_ES <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_ES, , file = here("Rdata", "MA_CVR_ES.rds"))Used to produce Fig. S5
MA_CVR_ES <- readRDS(here("Rdata", "MA_CVR_ES.rds"))
# telling ggplot to stop reordering factors
MA_CVR_ES$left_out <- as.factor(MA_CVR_ES$left_out)
MA_CVR_ES$left_out <- factor(MA_CVR_ES$left_out, levels = MA_CVR_ES$left_out)
# plotting
leaveoneout_ES <- ggplot(MA_CVR_ES) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_ES0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_ES0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_ES0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_ESdat$Study_ID <- as.integer(dat$Study_ID)Fig. S5 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the interaction of environmental enrichment x stress meta-analytic model using standardised mean difference (SMD) instead of lnRR
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -126.0571 252.1141 260.1141 270.1576 260.5793
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4653 0.6821 30 no Study_ID
## sigma^2.2 0.3698 0.6081 92 no ES_ID
## sigma^2.3 0.0000 0.0002 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 257.4673, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6880 0.1763 3.9017 91 0.0002 0.3377 1.0382 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 6.657130e-01 3.709364e-01 2.947766e-01 3.351907e-08
Risk of Bias
Including if the study was randomised as a moderator and conducting a pair-wise contrast
mod_randomised_ES <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~ROB_randomisation -
1, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_ES)
r2_ml(mod_randomised_ES)contra_mod_randomised_ES <- contrast_fun(data = dat, response = lnRR_ESa, moderator = ROB_randomisation,
VCV = VCV_ES)
res_table_mod_randomised_ES <- get_estimate(model = mod_randomised_ES, contra = contra_mod_randomised_ES,
moderator = ROB_randomisation)
res_table_mod_randomised_ES## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.0467 -0.106 0.200 0.546 -0.434 0.527
## 2 Yes 0.200 0.0503 0.350 0.00946 -0.279 0.680
## 3 Unclear-Yes 0.154 -0.0200 0.327 0.0821 -0.333 0.641
Including if the study used blinding and conducting a pair-wise contrast
mod_blinding_ES <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_ES)
r2_ml(mod_blinding_ES)contra_mod_blinding_ES <- contrast_fun(data = dat, response = lnRR_ESa, moderator = ROB_blinding,
VCV = VCV_ES)
res_table_mod_blinding_ES <- get_estimate(model = mod_blinding_ES, contra = contra_mod_blinding_ES,
moderator = ROB_blinding)
res_table_mod_blinding_ES## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.158 -0.000309 0.316 0.0504 -0.364 0.680
## 2 Yes 0.0444 -0.178 0.267 0.693 -0.500 0.589
## 3 Unclear-Yes -0.114 -0.343 0.115 0.327 -0.661 0.434
‘Pairwise’ effect sizes
We conducted five ‘meta-analyses on ’pair-wise’ comparisons using lnRR between the treatments and control: \[ \ln{\text{RR}} = \ln (M_{\text{EC}}/M_{\text{CC}}), \ln (M_{\text{CS}}/M_{\text{CC}}), \ln (M_{\text{ES}}/M_{\text{CC}}), \ln (M_{\text{ES}}/M_{\text{CS}}) \ln (M_{\text{ES}}/M_{\text{EC}}), \]
Enrichment relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{EC}}/M_{\text{CC}}) \]
VCV_E20 <- impute_covariance_matrix(vi = dat$lnRRV_E2, cluster = dat$Study_ID, r = 0.5)
mod_E20 <- rma.mv(yi = lnRR_E2a, V = VCV_E20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_E20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.3235 14.6470 22.6470 32.6904 23.1121
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0037 0.0611 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0281 0.1675 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 475.8327, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1066 0.0291 3.6655 91 0.0004 0.0489 0.1644 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.610789e-01 1.011724e-01 2.607945e-09 7.599065e-01
orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")Stress relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{CS}}/M_{\text{CC}}) \]
VCV_S20 <- impute_covariance_matrix(vi = dat$lnRRV_S2, cluster = dat$Study_ID, r = 0.5)
mod_S20 <- rma.mv(yi = lnRR_S2a, V = VCV_S20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_S20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -52.3561 104.7122 112.7122 122.7557 113.1773
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0323 0.1797 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0798 0.2824 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 1003.0694, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1846 0.0522 -3.5360 91 0.0006 -0.2882 -0.0809 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 9.489604e-01 2.734220e-01 9.879730e-10 6.755384e-01
orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{CC}}) \]
VCV_ES20 <- impute_covariance_matrix(vi = dat$lnRRV_ES2, cluster = dat$Study_ID,
r = 0.5)
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = VCV_ES20, random = list(~1 | Study_ID, ~1 |
Strain, ~1 | ES_ID), test = "t", data = dat)
summary(mod_ES20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.3625 20.7250 28.7250 38.7684 29.1901
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0039 0.0625 30 no Study_ID
## sigma^2.2 0.0014 0.0377 6 no Strain
## sigma^2.3 0.0227 0.1508 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 402.2656, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.0801 0.0389 2.0594 91 0.0423 0.0028 0.1573 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.81513970 0.11355557 0.04132363 0.66026050
orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to stress
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{CS}}) \]
VCV_E30 <- impute_covariance_matrix(vi = dat$lnRRV_E3, cluster = dat$Study_ID, r = 0.5)
mod_E30 <- rma.mv(yi = lnRR_E3a, V = VCV_E30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_E30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -46.3447 92.6895 100.6895 110.7329 101.1546
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0235 0.1532 30 no Study_ID
## sigma^2.2 0.0263 0.1623 6 no Strain
## sigma^2.3 0.0543 0.2331 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 790.0249, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.2465 0.1011 2.4389 91 0.0167 0.0457 0.4472 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.9456920 0.2132063 0.2391809 0.4933048
orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")Enrichment + stress relative to enrichment
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{EC}}) \]
VCV_S30 <- impute_covariance_matrix(vi = dat$lnRRV_S3, cluster = dat$Study_ID, r = 0.5)
mod_S30 <- rma.mv(yi = lnRR_S3a, V = VCV_S30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_S30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.8788 13.7576 21.7576 31.8011 22.2228
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0009 0.0292 6 no Strain
## sigma^2.3 0.0276 0.1662 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 540.3522, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0087 0.0342 -0.2552 91 0.7992 -0.0766 0.0592
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.415428e-01 5.192043e-09 2.515933e-02 8.163835e-01
orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")Figures
Panel of ‘focal’ ES and ‘pairwise’ ES orchard plots
Used to produce Fig. 3
mod_list1 <- list(mod_E0, mod_S0, mod_ES0)
mod_res1 <- lapply(mod_list1, function(x) mod_results(x, mod = "Int"))
merged1 <- submerge(mod_res1[[3]], mod_res1[[2]], mod_res1[[1]], mix = T)
merged1$mod_table$name <- factor(merged1$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
merged1$data$moderator <- factor(merged1$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
orchard1<- orchard_plot(merged1, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling +
scale_colour_manual(values = c("#00AEEF","#00A651","#ED1C24"))+ # change colours
scale_fill_manual(values=c("#00AEEF","#00A651","#ED1C24"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) mod_list2 <- list(mod_S30, mod_E30, mod_ES20, mod_S20, mod_E20) #rearranged the order so that it matches intext results
mod_res2 <- lapply(mod_list2, function(x) mod_results(x, mod = "Int"))
merged2 <- submerge(mod_res2[[1]], mod_res2[[2]], mod_res2[[3]], mod_res2[[4]], mod_res2[[5]], mix = T)
merged2$mod_table$name <- factor(merged2$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
merged2$data$moderator <- factor(merged2$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
orchard2 <- orchard_plot(merged2, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
scale_colour_manual(values = c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+ # change colours
scale_fill_manual(values=c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) p1 <- orchard1 + orchard2 + plot_annotation(tag_levels = "A")
p1# saved as PDF: 6 x 15 inchesFig. 3 Orchard plots showing the eight meta-analytic models. (A) main effects of enrichment and stress, and interaction of environmental enrichment and stress, (B) ‘pairwise’ comparisons between treatments and controls. Thick black lines = 95% CI, thin black lines = 95 % PI.
Panel of meta-regressions
Environmental enrichment
Used to produce Fig. 4
# Enrichment
E_mod <- (LvsM_E + Learning_E + Reinforcement_E)/(Age_E + Exercise_E + Social_E) +
plot_annotation(tag_levels = "A")
E_mod# saved as pdf 10 x 15 inchesFig. 4 Orchard plots showing the six different meta-regressions on the main effect of environmental enrichment on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at environmental enrichment, (E) type of manipulation of exercise during enrichment, (F) manipulation of the social environment during enrichment.
Stress
Used to produce Fig. 5
S_mod <- (LvsM_S + Learning_S + Reinforcement_S)/(Age_S + Stressor + Duration_S) +
plot_annotation(tag_levels = "A")
S_mod# saved as pdf 10 x 15 inchesFig. 5 Orchard plots showing the six different meta-regressions on the main effect of stress on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at stress, (E) the type of stressor (MS = maternal separation), (F) chronic or acute stress.
Interaction
Used to produce Fig. 6
ES_mod <- plot_grid(LvsM_ES, Learning_ES, Reinforcement_ES, Age_enrichment_ES, Age_stress_ES,
Order_ES, Exercise_ES, Social_ES, Stressor_ES, Duration_ES, labels = "AUTO",
ncol = 5)
ES_mod# saved as 10 x 25 inchesFig. 6 Orchard plots showing the 10 different meta-regressions of moderators on the interaction between environmental enrichment and stress learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement used, (D) the age at environmental enrichment, (E) the age at stress, (F) the order of treatment exposure, (G) if enrichment involved a manipulation of exercise, (H) manipulation of the social environment during enrichment, (I) the type of stressor, (F) stress was chronic or acute.
Panel of funnel plots
Used to produce Fig. 7
# EE
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
A <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
A <- recordPlot()
invisible(dev.off())
# Stress
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
B <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
B <- recordPlot()
invisible(dev.off())
# Interaction
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
C <- funnel(mod_ESfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
C <- recordPlot()
invisible(dev.off())
# putting together
ggdraw(A) + ggdraw(B) + ggdraw(C) + plot_annotation(tag_levels = "A")
# png(file = here('figs', 'Fig7_Funnels.png'))
# dev.off()knitr::include_graphics(here("figs", "funnels.png"))Fig. 7 Funnel plots of the standard error and residuals (lnRR) from the full models. (A) environmental enrichment main effect, (B) stress main effect, (C) environmental enrichment x stress interaction.
Software and package versions
sessionInfo() %>%
pander()R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: formatR(v.1.11), pander(v.0.6.4), gridGraphics(v.0.5-1), png(v.0.1-7), cowplot(v.1.1.1), ggthemr(v.1.1.0), ggalluvial(v.0.12.3), visdat(v.0.5.3), ggsignif(v.0.6.2), networkD3(v.0.4), GoodmanKruskal(v.0.0.3), patchwork(v.1.1.1), MuMIn(v.1.43.17), orchaRd(v.0.0.0.9000), clubSandwich(v.0.5.3), metafor(v.3.0-2), Matrix(v.1.3-3), here(v.1.0.1), forcats(v.0.5.1), stringr(v.1.4.0), dplyr(v.1.0.6), purrr(v.0.3.4), readr(v.1.4.0), tidyr(v.1.1.3), tibble(v.3.1.2), ggplot2(v.3.3.3) and tidyverse(v.1.3.1)
loaded via a namespace (and not attached): nlme(v.3.1-152), fs(v.1.5.0), lubridate(v.1.7.10), httr(v.1.4.2), rprojroot(v.2.0.2), tools(v.4.1.0), backports(v.1.2.1), utf8(v.1.2.1), R6(v.2.5.0), vipor(v.0.4.5), DBI(v.1.1.1), colorspace(v.2.0-1), withr(v.2.4.2), tidyselect(v.1.1.1), compiler(v.4.1.0), cli(v.2.5.0), rvest(v.1.0.0), pacman(v.0.5.1), xml2(v.1.3.2), sandwich(v.3.0-1), labeling(v.0.4.2), bookdown(v.0.22), scales(v.1.1.1), digest(v.0.6.27), rmarkdown(v.2.8), pkgconfig(v.2.0.3), htmltools(v.0.5.1.1), highr(v.0.9), dbplyr(v.2.1.1), htmlwidgets(v.1.5.3), rlang(v.0.4.11), readxl(v.1.3.1), rstudioapi(v.0.13), farver(v.2.1.0), generics(v.0.1.0), zoo(v.1.8-9), jsonlite(v.1.7.2), magrittr(v.2.0.1), Rcpp(v.1.0.6), ggbeeswarm(v.0.6.0), munsell(v.0.5.0), fansi(v.0.4.2), lifecycle(v.1.0.0), stringi(v.1.6.2), yaml(v.2.2.1), mathjaxr(v.1.4-0), crayon(v.1.4.1), lattice(v.0.20-44), haven(v.2.4.1), hms(v.1.1.0), knitr(v.1.33), pillar(v.1.6.1), igraph(v.1.2.6), codetools(v.0.2-18), stats4(v.4.1.0), reprex(v.2.0.0), glue(v.1.4.2), evaluate(v.0.14), modelr(v.0.1.8), vctrs(v.0.3.8), rmdformats(v.1.0.2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), xfun(v.0.23), broom(v.0.7.6), beeswarm(v.0.4.0) and ellipsis(v.0.3.2)
Social enrichment
Did enrichment involve more conspecifics (note that we did not include studies that only provided social enrichment but no other form of abiotic enrichment)?
Pairwise comparisons